Multi-class  Graph  Mumford-Shah  Model  for 
Plume  Detection  using  the  MBO  scheme* 


Huiyi  Hu^,  Justin  Sunu^,  and  Andrea  L.  Bertozzi^ 

^  Department  of  Mathematics,  University  of  California,  Los  Angeles 
huiyihuQmath . ucla . edu  bertozzi@math . ucla . edu 
^  Institute  of  Mathematical  Science,  Claremont  Graduate  University 
justinsunuQgmail . com 


Abstract.  We  focus  on  the  multi-class  segmentation  problem  using  the 
piecewise  constant  Mumford-Shah  model  in  a  graph  setting.  After  for¬ 
mulating  a  graph  version  of  the  Mumford-Shah  energy,  we  propose  an 
efficient  algorithm  called  the  MBO  scheme  using  threshold  dynamics. 
Theoretical  analysis  is  developed  and  a  Lyapunov  functional  is  proven 
to  decrease  as  the  algorithm  proceeds.  Furthermore,  to  reduce  the  com¬ 
putational  cost  for  large  datasets,  we  incorporate  the  Nystrom  extension 
method  which  efficiently  approximates  eigenvectors  of  the  graph  Lapla- 
cian  based  on  a  small  portion  of  the  weight  matrix.  Finally,  we  imple¬ 
ment  the  proposed  method  on  the  problem  of  chemical  plume  detection 
in  hyper-spectral  video  data. 

Keywords:  graph,  segmentation,  Mumford-Shah,  total  variation,  MBO, 
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1  Introduction 


Multi-class  segmentation  has  been  studied  as  an  important  problem  for  many 
years  in  various  areas,  such  as  computer  science  and  machine  learning.  For  im¬ 
agery  data  in  particular,  the  Mumford-Shah  model  [18]  is  one  of  the  most  exten¬ 
sively  used  model  in  the  past  decade.  This  model  approximates  the  true  image 
by  an  optimal  piecewise  smooth  function  through  solving  a  energy  minimiza¬ 
tion  problem.  More  detailed  review  of  the  work  on  Mumford-Shah  model  can  be 
found  in  the  references  of  [4].  A  simplified  version  of  Mumford-Shah  is  the  piece- 
wise  constant  model  (also  known  as  the  “minimal  partition  problem”),  which  is 
widely  used  due  to  its  reduced  complexity  compared  to  the  original  one.  For  a 
given  contour  ^  which  segments  an  image  region  i?  into  fi  many  disjoint  sub- 
regions  Q  =  the  piecewise  constant  Mumford-Shah  energy  is  defined 

as: 


(liQ  ^r)  j 


(1) 
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where  uq  is  the  observed  image  data,  {cr}^^i  is  a  set  of  constant  values,  and 
1^1  denotes  the  length  of  the  contour  By  minimizing  the  energy  over 

^  and  one  obtains  an  optimal  function  which  is  constant  within  each 

sub-region  to  approximate  uq,  along  with  a  segmentation  given  by  the  optimal 

In  [5],  a  method  of  active  contours  without  edges  is  proposed  to  solve  for 
the  two-class  piecewise  constant  Mumford-Shah  model  (n  =  2),  using  a  level  set 
method  introduced  in  [19].  The  work  in  [5]  is  further  generalized  to  a  multi-class 
scenario  in  [24].  The  method  developed  in  [5,  24]  is  well  known  as  the  Chan-Vese 
model,  which  is  a  popular  and  representative  method  for  image  segmentation. 
The  Chan-Vese  method  has  been  widely  used  due  to  the  model’s  flexibility  and 
the  great  success  it  achieves  in  performance. 

In  this  work,  we  formulate  the  piecewise  constant  MS  problem  in  a  graph 
setting  instead  of  a  continuous  one,  and  propose  an  efficient  algorithm  to  solve 
it.  Recently  the  authors  of  [2]  introduced  a  binary  semi-supervised  segmentation 
method  based  on  minimizing  the  Ginzburg-Landau  functional  on  a  graph.  In¬ 
spired  by  [2] ,  a  collection  of  work  has  been  done  on  graph-based  high-dimensional 
data  clustering  problems  posed  as  energy  minimization  problems,  such  as  semi- 
supervised  methods  studied  in  [14,11]  and  an  unsupervised  network  clustering 
method  [13]  known  as  modularity  optimization.  These  methods  make  use  of 
graph  tools  [6]  and  efficient  graph  algorithms,  and  our  work  pursues  similar 
ideas.  Note  that  unlike  the  Chan-Vese  model  which  uses  log2{n)  many  level 
set  functions  and  binary  representations  to  denote  multiple  classes,  our  model 
uses  simplex  constrained  vectors  for  class  assignments  representation  (details 
explained  below). 

To  solve  the  multi-class  piecewise  constant  MS  variational  problem  in  the 
graph  setting,  we  propose  an  efficient  algorithm  using  threshold  dynamics.  This 
algorithm  is  a  variant  of  the  one  presented  in  the  work  of  Merriman,  Bence  and 
Osher  (MBO)  [16, 17],  which  was  introduced  to  approximate  the  motion  of  an 
interface  by  its  mean  curvature  in  a  continuous  space.  The  idea  of  the  MBO 
scheme  is  used  on  the  continuous  MS  model  [8,  21]  motivated  by  level  set  meth¬ 
ods.  The  authors  of  [11, 13, 14]  implement  variants  of  the  MBO  scheme  applied 
to  segmentation  problems  in  a  graph  setting.  Rigorous  proofs  of  convergence  of 
the  original  MBO  scheme  in  continuous  setting  can  be  found  in  [1,9]  for  the 
binary  case,  and  [7]  for  the  multi-class  case.  An  analogous  discussion  in  a  graph 
setting  is  given  in  [23] .  Inspired  by  the  work  of  [7,23],  we  develop  a  Lyapunov 
functional  for  our  proposed  variant  of  the  MBO  algorithm,  which  approximates 
the  graph  MS  energy.  Theoretical  analysis  is  given  to  prove  that  this  Lyapunov 
energy  decreases  at  each  iteration  of  our  algorithm,  until  it  converges  within 
finitely  many  steps. 

In  order  to  solve  for  each  iteration  of  the  MBO  scheme,  one  needs  to  compute 
the  weight  matrix  of  the  graph  as  well  as  the  eigenvectors  of  the  corresponding 
graph  Laplacian.  However,  the  computational  cost  can  become  prohibitive  for 
large  datasets.  To  reduce  the  numerical  expenses,  we  implement  the  Nystrom 
extension  method  [10]  to  approximately  compute  the  eigenvectors,  which  only 
requires  computing  a  small  portion  of  the  weigh  matrix.  Thus  the  proposed 
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algorithm  is  efficient  even  for  large  datasets,  such  as  the  hyper-spectral  video 
data  considered  in  this  paper. 

The  proposed  method  can  be  implemented  on  general  high-dimensional  data 
clustering  problems.  However,  in  this  work  the  numerical  experiment  is  focused 
on  the  detection  of  chemical  plumes  in  hyper-spectral  video  data.  Detecting 
harmful  gases  and  chemical  plumes  has  wide  applicability,  such  as  in  environ¬ 
mental  study,  defense  and  national  security.  However,  the  diffusive  nature  of 
plumes  poses  challenges  and  difficulties  for  the  problem.  One  popular  approach 
is  to  take  advantage  of  hyper-spectral  data,  which  provides  much  richer  sensing 
information  than  ordinary  visual  images.  The  hyper-spectral  images  used  in  this 
paper  were  taken  from  video  sequences  captured  by  long  wave  infrared  (LWIR) 
spectrometers  at  a  scene  where  a  collection  of  plume  clouds  is  released.  Over 
100  spectral  channels  at  each  pixel  of  the  scene  are  recorded,  where  each  chan¬ 
nel  corresponds  to  a  particular  frequency  in  the  spectrum  ranging  from  7,820 
nm  to  11,700  nm.  The  data  is  provided  by  the  Applied  Physics  Laboratory  at 
Johns  Hopkins  University,  (see  more  details  in  [3]).  Prior  analysis  of  this  dataset 
can  be  found  in  the  works  [12, 15,  20,  22].  The  authors  of  [15]  implement  a  semi- 
supervised  graph  model  using  a  similar  MBO  scheme.  In  the  present  paper,  each 
pixel  is  considered  as  a  node  in  a  graph,  upon  which  the  proposed  unsuper¬ 
vised  segmentation  algorithm  is  implemented.  Competitive  results  are  achieved 
as  demonstrated  below. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  introduces  the  graph 
formula  for  the  multi-class  piecewise  constant  Mumford-Shah  model  and  relevant 
notations.  In  Secion  3,  the  Mumford-Shah  MBO  scheme  is  presented  as  well  as 
the  theoretical  analysis  for  a  Lyapunov  functional  which  is  proven  to  decrease  as 
the  algorithm  proceeds;  techniques  such  as  Nystrom  method  are  also  introduced 
for  the  purpose  of  numerical  efficiency.  In  Section  4,  our  algorithm  is  tested  on 
the  hyper-spectral  video  data  for  plume  detection  problem.  The  results  are  then 
presented  and  discussed. 


2  Graph  Mumford-Shah  Model 

Consider  an  A^-node  weighted  graph  (G,  L^),  where  G  =  {ui,  n2, . . . ,  is  a 
node  set  and  E  =  {wij}fLi  an  edge  set.  Each  node  rii  corresponds  to  an  agent 
in  a  given  dataset,  (such  as  a  pixel  in  an  image).  The  quantity  Wij  represents 
the  similarity  between  a  pair  of  nodes  and  rij.  Let  W  =  [wij]  denote  the 
graph’s  N  X  N  weight  matrix^  and  in  this  work  we  assume  W  is  symmetric,  i.e. 
Wij  =  Wji.  In  the  case  of  hyper- spectral  data,  each  node  (pixel)  rii  is  associated 
with  a  d-dimensional  feature  vector  (spectral  channels).  Let  uq  :  G  denote 

the  raw  hyper-spectral  data,  where  uo{ni)  represents  the  d-dimensional  spectral 
channels  of  rii.  We  use  the  following  notation: 

—  The  matrix  L  :=  D  —  W  is  called  the  (un- normalized)  graph  Laplacian 
[6],  where  D  is  a  diagonal  matrix  with  the  i-th  entry  being 
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:  G  ^  M,  observe  that 


1  ^ 

(t>,Lv)  =  2  • 

^i=i 

—  Graph  function  spaces  for  /  =  (/i,  /2,  •  •  • ,  fh)  •  G  ^ 

K:=|/|  /:G^[0,l]^E^(n,)  =  l| 

which  is  a  compact  and  convex  set. 


(2) 


B:=  |/|  /:G^{0,ir,^/rK)  =  l|  eK. 

This  simplex  constrained  vector  value  taken  by  /  G  ®  indicates  class  assign¬ 
ment,  i.e.  if  fr{rii)  =  1  for  some  r,  then  rii  belongs  to  the  r-th  class.  Thus 
for  each  /  G  ®,  it  corresponds  to  a  partition  of  the  graph  G  with  at  most  h 
classes.  Let  (/,L/)  =  E”=i(/r, L/^). 

—  Total  Variation  (TV)  for  graph  G  is  given  as: 

1  ^ 

I/Itv  •=  2  ^ 

^i=i 

In  this  setting,  we  present  a  graph  version  of  the  multi-class  piecewise  con¬ 
stant  Mumford-Shah  energy  functional: 


{G}r=l)  •“  2  ^ G|I  ^  fr)  1 

r=l 


(4) 


where  {cr}r=i  C  ||i4o  ~  denotes  an  N  x  1  vector 

{\\uo{n,)-Crf,...,\\MnN)-Crff  , 

and  {\\uo  —  c^|p,/r)  =  /r(^i)ll'^o(^i)  ~  Note  that  when  Ui  and  Uj 

belong  to  different  classes,  we  have  |/(^i)  — /(^j)|  =  2,  which  leads  to  the 
coefficient  in  front  of  the  term  ^  \f\j^y 

To  see  the  connection  between  (4)  and  (1),  one  first  observes  that  fr  is  the 
characteristic  function  of  the  r-th  class,  and  thus  {\\uq  —  c^|p,  fr)  is  analogous  to 
the  term  {uq  —  Cr)^  in  (1).  Furthermore,  the  total  variation  of  the  character¬ 
istic  function  of  a  region  gives  the  length  of  its  boundary  contour,  and  therefore 
is  the  graph  analogy  of  |^|. 

In  order  to  find  a  segmentation  for  G,  we  propose  to  solve  the  following 
minimization  problem: 


min 


MS{f,{Cr}%i). 


(5) 
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The  resulting  minimizer  /  yields  a  partition  of  G. 

One  can  observe  that  the  optimal  solution  of  (5)  must  satisfy: 


Cr 


{Upjr) 

T,f=l  fr{ni)  ’ 


(6) 


if  the  r-th  class  is  non-empty. 

Note  that  for  the  minimization  problem  given  in  (5),  it  is  essentially  equiva¬ 
lent  to  the  K-means  method  when  A  goes  to  +oo.  When  A  0,  the  minimizer 
approaches  a  constant. 


3  Mumford-Shah  MBO  and  Lyapunov  functional 

The  authors  of  [16, 17]  introduced  an  efficient  algorithm  (known  as  the  MBO 
scheme)  to  approximate  the  motion  by  mean  curvature  of  an  interface  in  a  con¬ 
tinuous  space.  The  general  procedure  of  the  MBO  scheme  alternates  between 
solving  a  linear  heat  equation  and  thresholding.  One  interpretation  of  the  scheme 
is  that  it  replaces  the  non-linear  term  of  the  Allen- Cahn  equation  with  thresh¬ 
olding  [8].  In  this  section  we  propose  a  variant  of  the  original  MBO  scheme  to 
approximately  find  the  minimizer  of  the  energy  MS{f,  {cr}r=i)  presented  in  (4). 
Inspired  by  the  work  of  [7,  23],  we  write  out  a  Lyapunov  functional  W(/)  for  our 
algorithm  and  prove  that  it  decreases  at  each  iteration  of  the  MBO  scheme. 


3.1  Mumford-Shah  MBO  scheme 


We  first  introduce  a  “diffuse  operator”  where  L  is  the  graph  Lapla- 

cian  defined  above  and  r  is  a  time  step  size.  The  operator  F^  is  analogous  to 
the  diffuse  operator  of  the  heat  equation  in  PDE  (continuous  space).  It 

satisfies  the  following  properties. 

Proposition  1.  Firstly,  F^  is  strictly  positive  definite,  i.e.  {f,F^f)  >  0  for  any 
/  G  K,  /  7^  0.  Secondly,  F^  conserves  the  mass,  i.e.  (1,/A/)  =  (1,/).  At  last, 
the  quantity  ^(1  —  f,F^f)  approximates  ^  Iflj^y,  for  any  /  G  ®. 

Proof.  Taylor  expansion  gives 


3"^^  =  I-tL+  +  •  •  • 

2!  3! 


Suppose  V  is  an  eigenvector  of  L  associated  with  the  eigenvalue  f.  One  then  has 
F^v  =  e~^^v  ^  {v,Frv)  =  e~^^{v,v)  >  0.  Let  the  eigen-decomposition  (with 
respect  to  L)  for  a  non-zero  /  :  G  ^  M  to  be  /  =  where  is 

a  set  of  orthogonal  eigenvectors  of  L  (note  that  L  is  positive  definite) .  Because 
//  is  a  linear  operator,  one  therefore  has  (/,  ///)  =  4>i)  >  0. 

For  the  second  property,  LI  =  0  ^  (1,L^/)  =  0,  where  1  is  an  N-dimensional 
vector  with  one  at  each  entry.  Therefore,  the  Taylor  expansion  of  /A  gives 

(i,r./)  =  (i,/). 
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At  last,  r,  c.  /  -  tL  ^  ^ (1  -  /,  r,/)  ~  i  (1  -  /,  /)  -  i  (1,  L/)  +  1  (/,  L/) . 
Particularly  when  /  G  B,  we  have  ^(1  —  /, /)  =  ^(1,L/)  =  0  and  |(/,  L/)  = 
^  l/l^y.  Hence  ^(1  -  f^r^f)  approximates  ^  |/|^^  in  B. 

□ 

Note  that  the  operator  (/  +  also  satisfies  the  above  three  properties,  and 

can  serve  the  same  purpose  as  as  far  as  this  paper  concerns. 

The  proposed  Mumford-Shah  MBO  scheme  for  the  minimization  problem  (5) 
consists  of  alternating  between  the  following  three  steps: 

For  a  given  G  B  at  the  k-th  iteration  and  , 

1.  Compute 

/  =  -  tX  (||mo  -  Cl  ||^  ||mo  -  clf,  •  •  •  ,  ||«0  -  c\f)  ,  (7) 

2.  (Thresholding) 

=  er,  r  =  aigma.xjc{ni) 

for  all  i  G  {1,  2, . . . ,  A^},  where  is  the  r-th  standard  basis  in  i.e. 
=  1  and  G  B. 

3.  (Update  c) 

fc+l  ^  juoJr^^) 


3.2  A  Lyapunov  functional 

We  introduce  a  Lyapunov  functional  W  for  the  Mumford-Shah  MBO  scheme: 


Yrif) 


^{l-f,rrf)+xf2{h0-Cr\\\fr), 

r=l 


subject  to  Cr  = 


Jr) 

JJr)  ■ 


According  to  the  third  property  of  U-  in  Proposition  1,  energy  Yr{f)  approxi¬ 
mates  MS{f for  /  G  B  and  A  similar  functional  for  the 

graph  total  variation  is  shown  and  discussed  in  [23]. 

Pursuing  similar  ideas  as  in  [7,23],  we  present  the  following  analysis  which 
consequently  shows  that  the  Mumford-Shah  MBO  scheme  (with  time  step  r) 
decreases  and  converges  to  a  stationary  state  within  a  finite  number  of  iter¬ 
ations. 

First  define 


Gr{f,  c)  :=  ^(1  -  f,  rrf)  +  A  ^(||wo  -  c.||2,  U) .  (9) 


Proposition  2.  The  functional  Gr{-,c)  is  strictly  concave  on  K,  for  any  fixed 
{c4”=i  e 
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Proof.  Take  /, ^  G  K,  a  G  (0, 1).  We  have  (1  —  a)f  +  G  K,  because  K  is  a 
convex  set. 

Gr  ((1  -  a)f  +  ag,  c)  -  (1  -  a)Gr{f,  c)  -  aGr{g,  c) 

=  ^a(l -«)(/- S'))  >0.  (10) 

Equality  only  holds  when  f  =  g.  Therefore,  Gr(',c)  is  strictly  concave  on  K. 

□ 

Aside  from  the  concavity  of  we  observe  that  the  first  order  variation  of 
Gr(',  c)  is  given  as 

fjGAf,c)  =  ^(1  -  ‘^r^f)  +  A  (lino  -  Cl||^  lino  -  C2||2,  . . .)  . 

Note  that  since  f)  is  linear,  the  Step  2  (thresholding)  in  the 

Mumford-Shah  MBO  scheme  is  equivalent  to 

Theorem  1.  In  the  Mumford-Shah  MBO  scheme,  the  Lyapunov  functional 
at  the  {k  +  l)-th  iteration  is  no  greater  than  Y^{f^).  Equality  only  holds  when 
fk  =  fk~^^.  Therefore,  the  scheme  achieves  a  stationary  point  inM  within  a  finite 
number  of  iterations. 

Proof. 

:=  argmin^e^(TG,(y,c'=),/)  (11) 

^  ^  to  linearity)  and 

0>{LGr{f,c^),f+^-f)  (12) 

>  o’")  -  Grip,  (concavity) 

Gr{P^^ ,c^)  <  Gt{P,c^)  =  Yr{P).  Observe  that  is  the 

minimizer  of 

^  GriP+\c'^+^)  <  Gr{P+\c'^)  <  YriP). 

^  Yr{fk~^^)  <  Yr{fk).  Therefore  the  Lyapunov  functional  W  is  decreasing  on 
the  iterations  of  the  Mumford-Shah  MBO  scheme,  unless  =  fk.  Since  ®  is 
a  finite  set,  a  stationary  point  can  be  achieved  in  a  finite  number  of  iterations. 

□ 

Minimizing  the  Lyapunov  energy  P^-  is  an  approximation  of  the  minimiza¬ 
tion  problem  in  (5),  and  the  proposed  MBO  scheme  is  proven  to  decrease  P^-. 
Therefore,  we  expect  the  Mumford-Shah  MBO  scheme  to  approximately  solve 
(5).  In  Section  3.3  and  Section  3.4,  we  introduce  techniques  for  computing  the 
MBO  iterations  efficiently. 


Huiyi  Hu,  Justin  Sunu,  and  Andrea  L.  Bertozzi 


3.3  Eigen-space  approximation 

To  solve  for  (7)  in  Step  1  of  the  Mumford-Shah  MBO  scheme,  one  needs  to 
compute  the  operator  T,-,  which  can  be  difficult  especially  for  large  datasets.  For 
the  purpose  of  efficiency,  we  numerically  solve  for  (7)  by  using  a  small  number  of 
the  leading  eigenvectors  of  L  (which  correspond  to  the  smallest  eigenvalues),  and 
project  onto  the  eigen-space  spanned  from  the  eigenvectors.  By  approximating 
the  operator  L  with  the  leading  eigenvectors,  one  can  compute  (7)  efficiently.  We 
use  this  approximation  because  in  graph  clustering  methods,  researchers  have 
been  using  a  small  portion  of  the  leading  eigenvectors  of  a  graph  Laplacian  to 
extract  structural  information  of  the  graph. 

Let  {0rn}m=i  dcuotc  the  first  M  (orthogonal)  leading  eigenvectors  of  L,  and 
corresponding  eigenvalues.  Assume  =  X]m=i  ^  where 

is  a  1 X  n  vector,  with  the  r-th  entry  =  (/^,  0m)-  Thus  /  can  be  approximately 
computed  as: 

M 

/  =  ^  -  tX  (||uo  -  c’lf,  ||wo  -  4^,  •  •  • ,  \\uo  -  )  .  (13) 

m=l 

The  Mumford-Shah  MBO  algorithm  with  the  above  eigen-space  approxima¬ 
tion  is  summarized  in  Algorithm  1.  After  the  eigenvectors  are  obtained,  each 
iteration  of  the  MBO  scheme  is  of  time  complexity  0{N).  Empirically,  the  al¬ 
gorithm  converges  after  a  small  number  of  iterations.  Note  that  the  iterations 
stop  when  a  purity  score  between  the  partitions  from  two  consecutive  iterations 
is  greater  than  99.9%.  The  purity  score,  as  used  in  [13],  measures  how  “similar” 
two  partitions  are.  Intuitively,  it  can  be  viewed  as  the  fraction  of  nodes  of  one 
partition  that  have  been  assigned  to  the  correct  class  with  respect  to  the  other 
partition. 


Algorithm  1  Mumford-Shah  MBO  algorithm 

Input:  /°,  Uo,  {(0m,^m}m=i,  r,  X,  n,  k  =  0. 
while  (purity (/^/^+^)  <  99.9%)  do 

-  /  =  Em=i  -  rX  (||i.o  -  ci|^  ||i.o  -  C2|^  . . . ,  ||i.o  -  c^H"). 

-  —  Gr,  where  r  —  argmax^ /c(ni). 

-  /c  ^  A;  +  1. 

end  while 


3.4  Nystrom  method 

The  Nystrom  extension  [10]  is  a  matrix  completion  method  which  has  been 
used  to  efficiently  compute  a  small  portion  of  the  eigenvectors  of  the  graph 
Laplacian  for  segmentation  problems  [2, 14, 11].  In  our  proposed  scheme,  leading 
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eigenvectors  of  L  are  required,  which  can  require  massive  computational  time 
and  memory.  For  large  graphs  such  as  the  ones  induced  from  images,  the  explicit 
form  of  the  weight  matrix  W  and  therefore  L  is  difficult  to  obtain  {0{N^)  time 
complexity).  Hence,  we  expect  to  use  the  Nystrom  method  to  approximately 
compute  the  eigenvectors  for  our  algorithm. 

Basically,  the  Nystrom  method  randomly  samples  a  very  small  number  (M) 
of  rows  of  W.  Based  on  matrix  completion  and  properties  of  eigenvectors,  it  ap¬ 
proximately  obtains  M  eigenvectors  of  the  symmetric  normalized  graph  Lapla- 
cian  L5  =  I  — D“2  WD“2  without  computing  the  whole  weight  matrix.  Detailed 
descriptions  of  the  Nystrom  method  can  be  found  in  [2, 14]. 

Note  that  our  previous  analysis  only  applies  to  L  rather  than  L^,  and  the 
Nystrom  method  can  not  be  trivially  formularized  for  L.  Therefore  this  question 
remains  to  be  studied.  However,  the  normalized  Laplacian  has  many  similar 
features  compared  to  L,  and  it  has  been  used  in  place  of  L  in  many  segmentation 
problems.  In  the  numerical  results  shown  below,  the  eigenvectors  of  computed 
via  Nystrom  perform  well  empirically. 

One  can  also  implement  other  efficient  methods  to  compute  the  eigenvectors 
for  the  Mumford-Shah  MBO  algorithm. 

4  Numerical  Results 


(a)  2nd  Eigenvector 


(b)  3rd  Eigenvector 


(c)  4th  Eigenvector 


(d)  5th  Eigenvector 


Fig.  1:  The  leading  eigenvectors  of  the  normalized  graph  Laplacian  computed  via 
the  Nystrom  method. 


The  hyper-spectral  images  tested  in  this  work  are  taken  from  the  video 
recording  of  the  release  of  chemical  plumes  at  the  Dugway  Proving  Ground, 
captured  by  long  wave  infrared  (LWIR)  spectrometers.  The  data  is  provided  by 
the  Applied  Physics  Laboratory  at  Johns  Hopkins  University.  A  detailed  descrip¬ 
tion  of  this  dataset  can  be  found  in  [3] .  We  take  seven  frames  from  a  plume  video 
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sequence  in  which  each  frame  is  composed  of  128  x  320  pixels.  We  use  a  back¬ 
ground  frame  and  the  frames  numbered  72  through  77  containing  the  plume. 
Each  pixel  has  129  spectral  channels  corresponding  to  a  particular  frequency 
in  the  EM  spectrum  ranging  from  7,820  nm  to  11,700  nm.  Thus,  the  graph  we 
construct  from  these  seven  frames  is  of  size  7  x  128  x  320  with  each  node  rii  cor¬ 
responding  to  a  pixel  with  a  129-dimensional  spectral  signature  Vi.  The  metric 
for  computing  the  weight  matrix  is  given  as: 


—  1 


which  is  an  approximation  of  the  spectral  angle  cos{vi,Vj). 

The  goal  is  to  segment  the  image  and  identify  the  “plume  cloud”  from  the 
background  components  (sky,  mountain, grass),  without  any  ground  truth.  As 
described  in  the  previous  section,  M  =  100  eigenvectors  of  the  normalized  graph 
Laplacian  (L^)  are  computed  via  the  Nystrom  method.  The  computational  time 
using  Nystrom  is  less  than  a  minute  on  a  2.8GHz  machine  with  Intel  Core  2Duo. 
The  visualization  of  the  first  five  eigenvectors  (associated  with  the  smallest  eigen¬ 
values)  are  given  in  Eigure  1  for  the  first  four  frames,  (the  first  eigenvector  is 
not  shown  because  it  is  close  to  a  constant  vector). 


(a) 


(b) 


Eig.  2:  The  segmentation  results  obtained  by  the  Mumford-Shah  MBO  scheme, 
on  a  background  frame  plus  the  frames  72-77.  Shown  in  (a)  and  (b)  are  segmen¬ 
tation  outcomes  obtained  with  different  initializations.  The  visualization  of  the 
segmentations  only  includes  the  first  four  frames. 


We  implement  the  Mumford-Shah  MBO  scheme  using  the  eigenvectors  on 
this  seven  frames  of  plume  images,  with  r  =  0.15,  A  =  150  and  n  =  5.  The 
test  is  run  for  20  times  with  different  uniformly  random  initialization,  and  the 
segmentation  results  are  shown  in  Eigure  2.  Note  that  depending  on  the  initial¬ 
ization,  the  algorithm  can  converge  to  different  local  minimum,  which  is  common 
for  most  non-convex  variational  methods.  The  result  in  (a)  occurred  five  times 
among  the  20  runs,  and  (b)  for  twice.  The  outcomes  of  other  runs  merge  either 
the  upper  or  the  lower  part  of  the  plume  with  the  background.  The  segmenta¬ 
tion  outcome  shown  in  (a)  gives  higher  energy  than  that  in  (b).  Among  the  20 
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runs,  the  lowest  energy  is  achieved  by  a  segmentation  similar  to  (a),  but  with 
the  lower  part  of  the  plume  merged  with  the  background.  It  may  suggest  that 
the  global  minimum  of  the  proposed  energy  does  not  necessarily  give  a  desired 
segmentation. 

Notice  that  in  Figure  2  (b),  even  though  there  actually  exist  five  classes,  only 
four  major  classes  can  be  perceived,  while  the  other  one  contains  only  a  very 
small  amount  of  pixels.  By  allowing  n  =  5  instead  of  4,  it  helps  to  reduce  the 
influence  of  a  few  abnormal  pixels.  The  computational  time  for  each  iteration 
is  about  2-3  seconds  on  a  1.7GHz  machine  with  Intel  Core  i5.  The  number  of 
iterations  is  around  20-40. 

3 

2.5 
2 

1.5 
1 

0.5 

Fig.  3:  Energy  MS{f)  (blue,  solid  line)  and  Y^{f)  (red,  dash  line)  at  each  itera¬ 
tion  from  the  same  test  as  shown  in  Figure  2  (a). 


Figure  3  demonstrates  a  plot  of  the  MS{f)  and  Yr{f)  energies  at  each  itera¬ 
tion  from  the  same  test  as  the  one  shown  in  Figure  2  (a).  The  Lyapunov  energy 
T,.(/)  (red,  dash  line)  is  non-increasing,  as  proven  in  Theorem  1.  Note  that  all 
the  energies  are  computed  approximately  using  eigenvectors. 

As  a  comparison,  the  segmentation  results  using  K-means  and  spectral  clus¬ 
tering  are  shown  in  Figure  4.  The  K-means  method  is  performed  directly  on 
the  raw  image  data  (7  x  128  x  320  by  129).  As  shown  in  (a)  and  (b),  the  re¬ 
sults  obtained  by  K-means  fail  to  capture  the  plume;  the  segmentations  on  the 
background  are  also  very  fuzzy.  For  the  spectral  clustering  method,  a  4- way  (or 
5- way)  K-means  is  implemented  on  the  four  (or  five)  leading  eigenvectors  of  the 
normalized  graph  Laplacian  (computed  via  Nystrom).  As  shown  in  (c)  and  (d), 
the  resulting  segmentations  divide  the  sky  region  into  two  undesirable  compo¬ 
nents.  Unlike  the  segmentation  in  Figure  2  (a)  where  the  mountain  component 
(red,  the  third  in  the  background)  has  a  well  defined  outline,  the  spectral  clus¬ 
tering  results  do  not  provide  clear  boundaries.  Our  approach  performs  better 
than  other  unsupervised  clustering  results  on  this  dataset  [12,20]. 

Another  example  of  the  plume  data  is  show  in  Figure  5,  where  the  67th  to 
72nd  frames  (instead  of  the  72nd  to  77th)  are  taken  along  with  the  background 
frame  as  the  test  data.  The  test  is  run  20  times  using  different  uniformly  random 
initialization,  where  r  =  0.15,  A  =  150  and  n  =  5.  The  result  in  Figure  5  (a) 
occurred  11  times  among  the  20  runs,  and  (b)  for  5  times.  The  outcomes  from 
the  other  4  runs  segment  the  background  into  three  components  as  in  (a),  but 
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(a)  4- way  K-means 


(b)  5-way  K-means 


(c)  Spectral  Clustering  with  4-way  K-means 


(d)  Spectral  Clustering  with  5-way  K-means 


Fig.  4:  K-means  and  spectral  clustering  segmentation  results.  The  visualization 
of  the  segmentations  only  includes  the  first  four  frames. 


Fig.  5:  The  segmentation  results  obtained  by  the  Mumford-Shah  MBO  scheme, 
on  a  background  frame  plus  the  frames  67-72.  Shown  in  (a)  and  (b)  are  segmen¬ 
tation  outcomes  obtained  with  different  initializations. 
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merge  the  plume  with  the  center  component.  The  segmentation  result  shown  in 
(a)  gives  the  lowest  energy  among  all  the  outcomes.  The  visualization  includes 
all  seven  frames  since  the  plume  is  small  in  the  first  several  frames. 

5  Conclusion 

In  this  paper  we  present  a  graph  framework  for  the  multi-class  piecewise  con¬ 
stant  Mumford-Shah  model  using  a  simplex  constrained  representation.  Based 
on  the  graph  model,  we  propose  an  efficient  threshold  dynamics  algorithm,  the 
Mumford-Shah  MBO  scheme  for  solving  the  minimization  problem.  Theoreti¬ 
cal  analysis  is  developed  to  show  that  the  MBO  iteration  decreases  a  Lyapunov 
energy  that  approximates  the  MS  functional.  Furthermore,  in  order  to  reduce 
the  computational  cost  for  large  datasets,  we  incorporate  the  Nystrom  exten¬ 
sion  method  to  approximately  compute  a  small  portion  of  the  eigenvectors  of 
the  normalized  graph  Laplacian,  which  does  not  require  computing  the  whole 
weight  matrix  of  the  graph.  After  obtaining  the  eigenvectors,  each  iteration  of 
the  Mumford-Shah  MBO  scheme  is  of  time  complexity  0{n).  The  number  of 
iterations  for  convergence  is  small  empirically. 

The  proposed  method  can  be  applied  to  general  high-dimensional  data  seg¬ 
mentation  problems.  In  this  work  we  focus  on  the  segmentation  of  hyper-spectral 
video  data.  Numerical  experiments  are  performed  on  a  collection  of  hyper- 
spectral  images  taken  from  a  video  for  plume  detection;  using  our  proposed 
method,  competitive  results  are  achieved.  However,  there  are  still  open  questions 
to  be  answered.  For  example,  the  Nystrom  method  can  only  compute  eigenvec¬ 
tors  for  the  normalized  Laplacian,  while  the  theoretical  analysis  for  the  Lyapunov 
functional  only  applies  to  the  un-normalized  graph  Laplacian.  This  issue  remains 
to  be  studied.  Note  that  the  graph  constructed  in  this  paper  does  not  include 
the  spacial  information  of  the  pixels,  but  only  the  spectral  information.  One 
can  certainly  build  a  graph  incorporating  the  location  of  each  pixel  as  well,  to 
generate  a  non-local  means  graph  as  discussed  in  [2]. 
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